#! /usr/bin/python

""" This script does a spectral analysis of a single data set and returns the power spectral density before, during and after a stimulus is presented. In addition a spectrogram is provided """


import os
import sys
import numpy
import scipy
import pylab

from LFP import signal_utils
from projects.electrophysiology.view_data import file_utils

filename = sys.argv[1]

time_resolution = 0.005 # in seconds
frequency_resolution = 1 # in Hz

stim_start = 10 #in secs 
stim_end = 15.3

def extract_data(filename):
    return file_utils.read_wessel_text_file(filename) 

voltage, sampling_freq, pulse_index = extract_data(filename)

stim_start_index = stim_start*sampling_freq
stim_end_index = stim_end*sampling_freq

#get the prestimulus PSD
#pre_psd_xx, pre_psd_freqs = signal_utils.psd(voltage[0:stim_start_index],
#                                     sampling_freq, 
#                                     frequency_resolution, 
#                                     high_frequency_cutoff=200)
#pylab.plot(pre_psd_freqs,pre_psd_xx, label='Prestimulus Power')
#pylab.hold(True)

#get the during stimulus PSD
#d_psd_xx, d_psd_freqs = signal_utils.psd(voltage[stim_start_index:stim_end_index],
#                                     sampling_freq, 
#                                     frequency_resolution, 
#                                     high_frequency_cutoff=200)

#pylab.plot(d_psd_freqs,d_psd_xx,label='During Stimulus Power')

#get the post stimulus PSD
#post_psd_xx, post_psd_freqs = signal_utils.psd(voltage[stim_end_index:len(voltage)],
#                                     sampling_freq, 
#                                     frequency_resolution, 
#                                     high_frequency_cutoff=200)
#pylab.plot(post_psd_freqs,post_psd_xx,label='Post Stimulus Power')

#pylab.xlim(1,100)
#pylab.ylim(0, .0003)
#pylab.legend(loc=7)
#pylab.show()

#pylab.subplot(212)
axes = pylab.gca()
Pxx, freqs, bins, im = signal_utils.specgram(voltage, sampling_freq, 
                        time_resolution, frequency_resolution, axes=axes, 
                        high_frequency_cutoff = 600,logscale=True)
pylab.ylim(0,100)
#pylab.xlabel('Time (s)', fontsize=15)
#pylab.ylabel('Frequency (Hz)', fontsize=15)
pylab.show()


